MODULE md1_multipole_ints

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_init_multipole_ints,           &
             md1_free_multipole_ints,           &
             md1_build_and_write_sph_mpoles

   REAL(REALK), ALLOCATABLE, SAVE :: car_mpoles(:,:,:)
   REAL(REALK), ALLOCATABLE, SAVE :: sph_mpoles(:,:,:)

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_init_multipole_ints(basis,LMAX)

      USE md1_car_to_sph, ONLY: md1_init_car_to_sph

      IMPLICIT NONE
      TYPE(md1_basis), INTENT(IN) :: basis
      INTEGER(INTK),   INTENT(IN) :: LMAX

      INTEGER(INTK) :: nAOmax

      nAOmax = (basis%LM(basis%Lmax)%car_nlm)**2
      ALLOCATE(sph_mpoles(nAOmax, (2*LMAX+1), 0:LMAX))
      ALLOCATE(car_mpoles(nAOmax, (LMAX+1)*(LMAX+2)/2, 0:LMAX))

      CALL md1_init_car_to_sph(LMAX)

   END SUBROUTINE md1_init_multipole_ints

!-------------------------------------------------------------------------------

   SUBROUTINE md1_free_multipole_ints

      USE md1_car_to_sph, ONLY: md1_free_car_to_sph

      IMPLICIT NONE
      CALL md1_free_car_to_sph
      IF (ALLOCATED(car_mpoles)) DEALLOCATE(car_mpoles)
      IF (ALLOCATED(sph_mpoles)) DEALLOCATE(sph_mpoles)

   END SUBROUTINE md1_free_multipole_ints

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_and_write_sph_mpoles(basis,LMAX,sh_pairs,nmoms)

      USE md1_car_to_sph, ONLY: md1_car_to_sph_mpoles

      IMPLICIT NONE
      TYPE(md1_basis),    INTENT(IN)  :: basis
      INTEGER(INTK),      INTENT(IN)  :: LMAX
      TYPE(md1_sh_pairs), INTENT(IN)  :: sh_pairs(:)
      INTEGER(INTK),      INTENT(OUT) :: nmoms

      INTEGER(INTK) :: ij, np

      nmoms = 0
      DO ij = 1, SIZE(sh_pairs)
         ! Build Cartesian multipoles for Cartesian GTO components
         CALL md1_build_car_mpoles(basis,sh_pairs(ij),LMAX,np)
         ! Transform Cartesian multipoles to spherical multipoles
         CALL md1_car_to_sph_mpoles(car_mpoles(:np,:,0:),               &
                                    sph_mpoles(:np,:,0:),np,LMAX)
         ! Transform Cartesian GTO components to spherical components
!         CALL md1_car_to_sph_GTOs(car_mpoles(:np,:,0:),                 &
!                                  sph_mpoles(:np,:,0:),np,LMAX)
         ! Write fully spherical integrals and overlap matrix to disk
         CALL md1_write_sph_mpoles(sh_pairs(ij),LMAX,nmoms)
      END DO

   END SUBROUTINE md1_build_and_write_sph_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_car_mpoles(basis,sh,LMAX,np)

      USE md1_utils,          ONLY: md1_get_prim_batch
      USE md1_hermite_mpoles, ONLY: md1_build_1D_mpoles

      IMPLICIT NONE
      TYPE(md1_basis),    INTENT(IN)  :: basis
      TYPE(md1_sh_pairs), INTENT(IN)  :: sh
      INTEGER(INTK),      INTENT(IN)  :: LMAX
      INTEGER(INTK),      INTENT(OUT) :: np

      REAL(REALK), ALLOCATABLE :: Qjie(:,:,:,:,:)
      TYPE(md1_prim_batch), POINTER :: batch(:)
      TYPE(md1_prim_batch), TARGET  :: tmp(basis%Cmax2)
      INTEGER(INTK) :: xyz

      ! Prepare batch of primitives, including multipole expansion origin
      CALL md1_get_prim_batch(basis,sh,tmp,np)
      batch => tmp(1:np)

      ! Build 1-D multipole integrals
      ALLOCATE (Qjie(np, 0:sh%Jang, 0:sh%Iang, 0:LMAX, 3))
      DO xyz = 1, 3
         CALL md1_build_1D_mpoles(xyz,batch,sh%xyz,Qjie(:,:,:,:,xyz))
      END DO
!print *, "--------------"
!print *, "shell pair centre:", sh%xyz
!do xyz=1, np
!   print *, xyz
!   print *, batch(xyz)%P
!   print *, batch(xyz)%PA
!   print *, batch(xyz)%PB
!   print *, batch(xyz)%exp
!   print *, batch(xyz)%KAB
!   print *, batch(xyz)%CAB
!end do

      ! Multiply 1-D integrals to form 3-D cartesian multipoles
      CALL md1_build_3D_mpoles(Qjie,basis%LM,batch,sh%I,sh%J)
      np = basis%LM(sh%Iang)%car_nlm * basis%LM(sh%Jang)%car_nlm
      DEALLOCATE (Qjie)

   END SUBROUTINE md1_build_car_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_3D_mpoles(Qjie,LM,batch,I,J)

      IMPLICIT NONE
      REAL(REALK),  TARGET, INTENT(IN) :: Qjie(:,0:,0:,0:,:)
      TYPE(md1_basis_LM),   INTENT(IN) :: LM(0:)
      TYPE(md1_prim_batch), INTENT(IN) :: batch(:)
      INTEGER(INTK),        INTENT(IN) :: I,J

      REAL(REALK), POINTER :: X(:), Y(:), Z(:)
      INTEGER(INTK) :: ituv(3), jtuv(3), etuv(3)
      INTEGER(INTK) :: Iang, Jang, L, pM, qM, qMhi, M, PQ

      Iang = UBOUND(Qjie,3)
      Jang = UBOUND(Qjie,2)

      car_mpoles(:,:,:) = zero

!print *, "car_mpoles: (L,M,PQ)"
      DO L = 0, UBOUND(Qjie,4)
         DO M = 1, LM(L)%car_nlm
            etuv(:) = LM(L)%tuv(:,M)

            PQ = 0
            DO pM = 1, LM(Iang)%car_nlm
               ituv(:) = LM(Iang)%tuv(:,pM)
               qMhi = LM(Jang)%car_nlm
               IF (I == J) qMhi = pM
               DO qM = 1, qMhi
                  PQ = PQ + 1
                  jtuv(:) = LM(Jang)%tuv(:,qM)
   
                  X => Qjie(:,jtuv(1),ituv(1),etuv(1),1)
                  Y => Qjie(:,jtuv(2),ituv(2),etuv(2),2)
                  Z => Qjie(:,jtuv(3),ituv(3),etuv(3),3)
                  car_mpoles(PQ,M,L) = SUM( X(:)*Y(:)*Z(:) * batch(:)%CAB )
!print *, L,M,PQ,car_mpoles(PQ,M,L)

               END DO
            END DO

         END DO
      END DO

   END SUBROUTINE md1_build_3D_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_write_sph_mpoles(sh,LMAX,nmoms)

      IMPLICIT NONE
      TYPE(md1_sh_pairs), INTENT(IN)    :: sh
      INTEGER(INTK),      INTENT(IN)    :: LMAX
      INTEGER(INTK),      INTENT(INOUT) :: nmoms

      REAL(REALK), PARAMETER :: MomScrn = 1d-15
      INTEGER(INTK) :: L, M, PQ, P,Q, Qhi
      LOGICAL :: flag

      PQ = 0
      DO P = sh%P%lo, sh%P%hi
         Qhi = sh%Q%hi
         IF (sh%I == sh%J) Qhi = P
         DO Q = sh%Q%lo, Qhi
            PQ = PQ + 1

            flag = .FALSE.
            DO L = 0, LMAX
               DO M = 1, 2*L + 1
                  ! Weak screening
                  IF (ABS(sph_mpoles(PQ,M,L)) < MomScrn) CYCLE
                  flag = .TRUE.
                  ! Write-out overlap matrix for debugging (L=0)
                  IF (L == 0) WRITE(IO3,'(3I8,F25.15)')             &
                              nmoms+1, P,Q, sph_mpoles(PQ,M,L)
                  ! Write multipoles to file
                  WRITE(IO2) nmoms+1, L, (M-L-1), P, Q,            &
                             sh%xyz(1:3), sph_mpoles(PQ,M,L)
                  WRITE(IO1,'(5I4,4F12.5)') nmoms+1, L, (M-L-1), P, Q,    &
                                 sh%xyz(1:3), sph_mpoles(PQ,M,L)
               END DO
            END DO
            IF (flag) nmoms = nmoms + 1

         END DO
      END DO

   END SUBROUTINE md1_write_sph_mpoles

!-------------------------------------------------------------------------------

END MODULE md1_multipole_ints

